function [rho_save,j,s,dWk,Am, int_field_000, int_field_011] = get_trajectory(init,T,dt, gamma_1,  gamma_d1, gamma_d2, gamma_d3, dim , kappa1, kappa2, epsilon0, chi, g, Delta_q, N, ton, toff, sigma, phi, tmeas_on, tmeas_off, Delta_a, Delta_b,eta )

%global gamma_d1 gamma_d2 gamma_d3 dim kappa1 kappa2 epsilon0 chi g Delta_q ton sigma phi tmeas_on


N = round(T/dt);                  % number of steps
s = 0;
dWk = zeros(1,N);

% Allocate memory to store the data
rho_save = zeros(dim^2,N);           % the density matrix
j = zeros(1,N-1);                    % the current

% Initial state
rho = init;
rho_save(:,1) = init;

sigmaZ = sparse([1,0;0,-1]);
sigmaM = sparse([0,1;0,0]);

eye = sparse(diag(ones(2,1)));
fulleye = sparse(diag(ones(dim,1)));


% The dephasing operator for qubit 1

%Ad1 = mult(sigmaZ,eye,eye); 
ZII = [ones(1,4),-ones(1,4)];
%Ad1 = diag(ZII);

%Ld1 = kron(Ad1,transpose(Ad1)) - 0.5*kron(Ad1'*Ad1,fulleye) - 0.5*kron(fulleye,transpose(Ad1'*Ad1));
ZII4 = [ZII, ZII, ZII, ZII];
Ld1 = spdiags([ZII4'; -ZII4']-1, 0, 64, 64); 

% The dephasing operator for qubit 2

%Ad2 = mult(eye,sigmaZ,eye);
IZI = [1, 1, -1, -1, 1, 1, -1, -1];
IZI2 = [IZI,IZI];
IZI8 = [IZI2, -IZI2, IZI2, -IZI2];

%Ld2 = kron(Ad2,transpose(Ad2)) - 0.5*kron(Ad2'*Ad2,fulleye) - 0.5*kron(fulleye,transpose(Ad2'*Ad2));
Ld2 = spdiags(IZI8'-1, 0, 64, 64);

% The dephasing operator for qubit 3

%Ad3 = mult(eye,eye,sigmaZ);
IIZ = [1, -1, 1, -1, 1, -1, 1, -1];
IIZ2 = [IIZ, -IIZ];
IIZ8 = [IIZ2, IIZ2, IIZ2, IIZ2];

%Ld3 = kron(Ad3,transpose(Ad3)) - 0.5*kron(Ad3'*Ad3,fulleye) - 0.5*kron(fulleye,transpose(Ad3'*Ad3));
Ld3 = spdiags(IZI8'-1, 0, 64, 64);

% The relaxation operator for qubit 1

%Ar1 = mult(sigmaM,eye,eye);
Ar1 = sparse(1:4,5:8,ones(4,1),8,8);
I4Ar1 = sparse(blkdiag(Ar1,Ar1,Ar1,Ar1));
SminusI4Ar1 = sparse([zeros(32), I4Ar1; zeros(32), zeros(32)]); %FIXME Rename this something legible.
%Ar1dagger = Ar1';
alt_zero_1_4_4 = [zeros(4,1);ones(4,1)];
I8kronAr1ish = 0.5*[alt_zero_1_4_4; alt_zero_1_4_4; alt_zero_1_4_4; alt_zero_1_4_4]; %FIXME this name is awful.
Lr1 = SminusI4Ar1 - spdiags([zeros(32,1);0.5*ones(32,1)],0, 64, 64) - spdiags(I8kronAr1ish, 0, 64, 64);

% The relaxation operator for qubit 2

sigmaM_I = [zeros(2) eye(2); zeros(2) zeros(2)];
I4_sigmaM_I = blkdiag(sigmaM_I,sigmaM_I,sigmaM_I,sigmaM_I);
sigmaM_I4_sigmaM_I = sparse([zeros(32) I4_sigmaM_I; zeros(32) zeros(32)]);

s1_I_d = [0;0;0.5;0.5];
I4_s1_I_d = [s1_I_d;s1_I_d;s1_I_d;s1_I_d];
I16_s1_I_d = [I4_s1_I_d;I4_s1_I_d;I4_s1_I_d;I4_s1_I_d];
%Ar2dagger = Ar2';
%Lr2 = kron(Ar2,transpose(Ar2dagger)) - 0.5*kron(Ar2dagger*Ar2,fulleye) - 0.5*kron(fulleye,transpose(Ar2dagger*Ar2));
Lr2 = sparse(blkdiag(sigmaM_I4_sigmaM_I,sigmaM_I4_sigmaM_I)) - spdiags([zeros(16,1),0.5*ones(16,1),zeros(16,1),0.5*ones(16,1)],0,64,64)...
      - spdiags(I16_s1_I_d,0,64,64);

% The relaxation operator for qubit 3

%Ar3 = mult(eye,eye,sigmaM);
Ar3 = spdiags([0,0,0,0,1,1,1,1],0,8,8);
%Ar3dagger = Ar3';
sm_Ar3 = sparse([zeros(8) Ar3; zeros(8) zeros(8)]);
alt8 = [zeros(8,1);0.5*ones(8,1)];
alt2 = [0,0.5,0,0.5,0,0.5,0,0.5,0,0.5,0,0.5,0,0.5,0,0.5]';
Lr3 = sparse(blkdiag(sm_Ar3,sm_Ar3,sm_Ar3,sm_Ar3)) - spdiags([alt8,alt8,alt8,alt8],0,64,64)...
    - spdiags([alt2;alt2;alt2;alt2],0,64,64);

% The time independent deterministic part of Liouvillian

D1 = 0.5*gamma_d1*Ld1 + 0.5*gamma_d2*Ld2 + 0.5*gamma_d3*Ld3 + 0*((g/Delta_q)^2*(kappa1+kappa2) + gamma_1 )*(Lr1 + Lr2 + Lr3);


%% Calculate the chi-matrix elements

chi_mat = zeros(2,2,2);

for a = [0,1]
    
    for b = [0,1]
    
        for c = [0,1]
            
            chi_mat(a+1,b+1,c+1) = -chi*((-1)^a + (-1)^b + (-1)^c);

        end
        
    end
    
end


%% Calculate the fields in the cavity

tvec = linspace(0,T,N);
%epsilon = epsilon0/pi*( atan(sigma*(tvec - ton)) + 0.5*pi);

epsilon = epsilon0/pi*( atan(sigma*(tvec - ton)) + 0.5*pi) - ...
          epsilon0/pi*( atan(sigma*(tvec - toff)) + 0.5*pi);

alpha = zeros(N,2,2,2);
beta = zeros(N,2,2,2);

for a = [0,1]
    
    for b = [0,1]
    
        for c = [0,1]
            
            [~,sol] = ode45(@(t,y) fieldeqns(t, y, tvec, epsilon, chi, kappa1, kappa2, Delta_a, Delta_b, a, b, c), tvec, [0,0]');
            
            alpha(:,a+1,b+1,c+1) = sol(:,1);
            beta(:,a+1,b+1,c+1)  = sol(:,2);
            
        end
        
    end
    
end

int_field_000 = trapz(tvec, sqrt(kappa1)*(alpha(:,1,1,1) + beta(:,1,1,1)) ); 

int_field_011 = trapz(tvec, sqrt(kappa1)*(alpha(:,1,2,2) + beta(:,1,2,2)) );

%% Calculate the projection operators 

Pi = zeros(dim,dim, 2,2,2);

for a = [0,1]
    
    for b = [0,1]
    
        for c = [0,1]
            
            indx = 4*a + 2*b + 1*c + 1;
            Pi(indx,indx,a+1,b+1,c+1) = 1;
            Pi(:,:,a+1,b+1,c+1) = sparse(Pi(:,:,a+1,b+1,c+1));
            
            
        end
        
    end
    
end



%% Time stepping using Kloeden-Platen method with fixed stepsize dt

for i=1:N-1
    
    % Generate a random number
    dWk(i) = sqrt(dt)*randn(1);
        
    % the time dependent measurement operator
    
    Am = zeros(dim,dim);
    
    for a = [0,1]
        
        for b = [0,1]
            
            for c = [0,1]
                
                Am = Am + ( sqrt(kappa1*eta)*alpha(i,a+1,b+1,c+1) + sqrt(kappa2*eta)*beta(i,a+1,b+1,c+1) ) * Pi(:,:,a+1,b+1,c+1);
                
            end
            
        end
        
    end
    
    Am = ( heaviside(i*dt - tmeas_on) - heaviside(i*dt - tmeas_off) )*exp(1j*phi)*Am;
    
    D2 = kron(Am,fulleye) + kron(fulleye,transpose(conj(Am)));
    
    % the current
    
    average = trace( ( Am + conj(Am) )*transpose(reshape(rho,dim,dim)));
    j(i) = average + dWk(i)/dt;
    
    % the integrated current
    s = s + average * dt + dWk(i);
    
    % The time dependent part of the deterministic super operator
    
    D1_t = zeros(dim^2,dim^2);
    
    for a = [0,1]
        
        for b = [0,1]
            
            for c = [0,1]
                
                for d = [0,1]
                    
                    for e = [0,1]
                    
                        for f = [0,1]
                
                            D1_t = D1_t +   ( (chi_mat(d+1,e+1,f+1) - chi_mat(a+1,b+1,c+1)) * imag(conj(alpha(i, a+1, b+1, c+1))*alpha(i, d+1, e+1, f+1) ) + ...
                                              (chi_mat(d+1,e+1,f+1) - chi_mat(a+1,b+1,c+1)) * imag(conj( beta(i, a+1, b+1, c+1))* beta(i, d+1, e+1, f+1)  ) ) ... 
                                              *kron(Pi(:,:, a+1, b+1, c+1), Pi(:,:, d+1, e+1, f+1));
                            
                        end
                        
                    end
                    
                end
                
            end
            
        end
        
    end

    for a = [0,1]
        
        for b = [0,1]
            
            for c = [0,1]
                            
                d = mod(a + 1,2);
                e = mod(b + 1,2);
                f = mod(c + 1,2);
                
                D1_t = D1_t + sqrt(kappa1*kappa2)*...
                                            ( real( beta(i, a+1, b+1, c+1)*conj(alpha(i, d+1, e+1, f+1) )) + ...
                                              real(alpha(i, a+1, b+1, c+1)*conj( beta(i, d+1, e+1, f+1) )) - ...
                                              real( beta(i, a+1, b+1, c+1)*conj(alpha(i, a+1, b+1, c+1) )) - ...
                                              real( beta(i, d+1, e+1, f+1)*conj(alpha(i, d+1, e+1, f+1) )) ) ... 
                                              *kron(Pi(:,:, a+1, b+1, c+1), Pi(:,:, d+1, e+1, f+1));
                
            end
            
        end
        
    end
    
    %% Order 2 weak Runge-Kutta scheme (Platen)
    
    
    d1 = (D1 + D1_t)*rho;                % D1(psi)
    
    
    d2 = (D2 - kron(fulleye,fulleye)*average)*rho; % D2(psi)
    
    zt=rho+d1*dt+d2*dWk(i);      % psi_tilde
    zp=rho+d1*dt+d2*sqrt(dt); % psi_plus
    zm=rho+d1*dt-d2*sqrt(dt); % psi_minus
    
    
    d1t = (D1 + D1_t)*zt;                % D1(psi)
    
    
    average_p = trace( (Am + conj(Am) )*transpose(reshape(zp,dim,dim)));
    average_m = trace( (Am + conj(Am) )*transpose(reshape(zm,dim,dim)));
    
    d2p = (D2 - kron(fulleye,fulleye)*average_p)*zp;
    d2m = (D2 - kron(fulleye,fulleye)*average_m)*zm;
    
    rho = rho + 0.5 * ( d1t + d1 ) * dt ...
        + 0.25 * ( d2p + d2m + 2*d2 ) * dWk(i) ...
        + 0.25 * ( d2p - d2m ) * ( dWk(i)^2 - dt ) / sqrt(dt);
    
    %rho_save_anal(:,i+1) = expm(0.5*phase_det*D)*expm(0*phase_stoch(i+1)*D2)*init;
    rho_save(:,i+1) = rho;
    
    
    
    %% Order 1 strong scheme (Milstein)
    %
    %         d1 = (D1 + D1_t)*rho;
    %         d2 = (D2 - sqrt(kappa)*kron(fulleye,fulleye)*average)*rho;
    %
    %         G = rho + d2*sqrt(dt);
    %         average_G = trace( (Am + conj(Am) )*reshape(G,4,4));
    %         d2G = (D2 - sqrt(kappa)*kron(fulleye,fulleye)*average_G)*G;
    %
    %         rho = rho + 0*d1*dt + 1*(d2*dWk + 0.5*( d2G - d2 )*( dWk^2 - dt ) / sqrt(dt));
    %         rho_save(:,i+1) = rho;
    
end
%s_phase = angle(rhoe(7));
% subtract the accumulated phase
%rho_save(:,end) = expm(0.5*s_phase*D_undo)*rho_save(:,end);

end



